Stroke vs. Bleeding in MIDAS Patients within a Year of Initial AF Diagnosis

1 Stroke vs. Bleeding Shiny App Design Proposal

# To Do 1. Work with a stats student (Yayie Duan) to write a Shiny app.

2. After meeting with Javier (01/28/2021):
   

1. Instead of linear regression, use Deming regression
   
2. Probability of being in the green region based on the confidence intervals (taken care of Deming regression?)
   
3. Think about: ratio of variance (related to Deming regression)

3. After meeting with Cande (01/28/2021):
   

1. Add more variables to the model. I’m thinking of adding CHADSVASC variables and displaying the score additional to the recommendation. We can also extend it to a more comprehensive model and use the app to also demonstrate how the MIDAS-based risk calculator compares to CHADSVASC risk estimation.
   b. IF we are going to assign probabilities based on confidence intervals (see ‘b’ above), we can have an internal rule, e.g. if p > 0.5 then medicate, if < then don’t). This does not need to be displayed in the app to make it less confusing for the patients.

2 Setup

require(data.table)

Loading required package: data.table
Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     
data.table 1.13.6 using 2 threads (see ?getDTthreads).  Latest news: r-datatable.com

require(ggplot2)

Loading required package: ggplot2

require(plotly)

Loading required package: plotly
Registered S3 methods overwritten by 'htmltools':
  method               from         
  print.html           tools:rstudio
  print.shiny.tag      tools:rstudio
  print.shiny.tag.list tools:rstudio
Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio

Attaching package: 㤼㸱plotly㤼㸲

The following object is masked from 㤼㸱package:ggplot2㤼㸲:

    last_plot

The following object is masked from 㤼㸱package:stats㤼㸲:

    filter

The following object is masked from 㤼㸱package:graphics㤼㸲:

    layout

require(deming)

Loading required package: deming

require(MethComp)

Loading required package: MethComp

ATTENTION! In ggplot, must explicitly set ‘inherit.aes = FALSE’, otherwise ‘geom_ribbon’ expects values for aesthetics from previous geoms (i.e. tt variables ps, pb, etc.). See solution here:
(geom_ribbon expects an unused y-axis variable to be present)[https://github.com/tidyverse/ggplot2/issues/3364]

3 Load data

load("data/tt.RData")

4 Metod 1: Linear regression

4.1 Run ‘lm’

# Log2-transformed variables
tt[, c("l2pb",
       "l2ps") := list(log2(pb), #xx
                        log2(ps))] #yy
# Regression on log2-transformed variables
m2 <- lm(l2ps ~ l2pb,
         data = tt)
summary(m2)

Call:
lm(formula = l2ps ~ l2pb, data = tt)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.67656 -0.50577  0.02698  0.40965  0.72161 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -4.5353     1.4238  -3.185  0.00973 ** 
l2pb          1.4596     0.1978   7.379 2.37e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.553 on 10 degrees of freedom
Multiple R-squared:  0.8449,    Adjusted R-squared:  0.8293 
F-statistic: 54.45 on 1 and 10 DF,  p-value: 2.371e-05

new_l2pb <- seq(from = 5.5,
                to = 8.75,
                length.out = nrow(tt))
# new_l2pb <- seq(from = 5, 
#                 to = 9,
#                 length.out = nrow(tt))
pi_20_80 <- data.table(predict(m2, 
                               newdata = list(l2pb = new_l2pb),
                               interval = "prediction",
                               level = 0.6))
pi_20_80$l2pb <- new_l2pb
head(pi_20_80)

4.2 Plot Fear Bleeding 20 on 0-100 Scale (Log2)

p1 <- ggplot(tt,
             aes(x = l2pb,
                 y = l2ps,
                 fill = lab)) +
  geom_errorbar(aes(ymin = log2(10^4*(prob.stroke - 1.96*se.stroke)),
                    ymax = log2(10^4*(prob.stroke + 1.96*se.stroke))),
                size = 0.5) +
  geom_errorbar(aes(xmin = log2(10^4*(prob.bleed - 1.96*se.bleed)),
                    xmax = log2(10^4*(prob.bleed + 1.96*se.bleed))),
                size = 0.5) +
  geom_point(shape = 21,
             aes(size = sqrt(N))) +
  geom_abline(slope = m2$coefficients[2],
              intercept = m2$coefficients[1],
              linetype = "dashed") +
  geom_ribbon(data = pi_20_80,
              aes(x = l2pb,
                  ymin = lwr,
                  ymax = 9.25),
              fill = "green",
              alpha = 0.2,
              inherit.aes = FALSE) + 
  scale_x_continuous("Log2 of Bleeding per 10^4 patients",
                      # limits = c(5.5, 8.5),
                     expand = c(0, 0)) +
  scale_y_continuous("Log2 of Stroke per 10^4 patients",
                     # limits = c(3, 9),
                     expand = c(0, 0)) +
  ggtitle("Fear of Bleeding: 20/100\nGREEN = MEDICATE") +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5))
ggplotly(p1)

4.3 Plot Fear Bleeding 80 on 0-100 Scale (Log2)

p2 <- ggplot(tt,
             aes(x = l2pb,
                 y = l2ps,
                 fill = lab)) +
  geom_errorbar(aes(ymin = log2(10^4*(prob.stroke - 1.96*se.stroke)),
                    ymax = log2(10^4*(prob.stroke + 1.96*se.stroke))),
                size = 0.5) +
  geom_errorbar(aes(xmin = log2(10^4*(prob.bleed - 1.96*se.bleed)),
                    xmax = log2(10^4*(prob.bleed + 1.96*se.bleed))),
                size = 0.5) +
  geom_point(shape = 21,
             aes(size = sqrt(N))) +
  geom_abline(slope = m2$coefficients[2],
              intercept = m2$coefficients[1],
              linetype = "dashed") +
  geom_ribbon(data = pi_20_80,
              aes(x = l2pb,
                  ymin = upr,
                  ymax = 9.25),
              fill = "green",
              alpha = 0.2,
              inherit.aes = FALSE) + 
  scale_x_continuous("Log2 of Bleeding per 10^4 patients",
                     expand = c(0, 0)) +
  scale_y_continuous("Log2 of Stroke per 10^4 patients",
                     expand = c(0, 0)) +
  ggtitle("Fear of Bleeding: 80/100\nGREEN = MEDICATE") +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5))
ggplotly(p2)

4.4 Readjust predictions

new_l2pb <- seq(from = 5,
                to = 9,
                length.out = 10*nrow(tt))
pi_20_80 <- data.table(predict(m2, 
                               newdata = list(l2pb = new_l2pb),
                               interval = "prediction",
                               level = 0.6))
pi_20_80$l2pb <- new_l2pb
head(pi_20_80)

4.5 Plot Fear Bleeding 20 on 0-100 Original Scale

p1 <- ggplot(tt,
             aes(x = pb,
                 y = ps,
                 fill = lab)) +
  geom_errorbar(aes(ymin = 10^4*(prob.stroke - 1.96*se.stroke),
                    ymax = 10^4*(prob.stroke + 1.96*se.stroke)),
                size = 0.5) +
  geom_errorbar(aes(xmin = 10^4*(prob.bleed - 1.96*se.bleed),
                    xmax = 10^4*(prob.bleed + 1.96*se.bleed)),
                size = 0.5) +
  geom_point(shape = 21,
             aes(size = sqrt(N))) +
  geom_ribbon(data = pi_20_80,
              aes(x = 2^l2pb,
                  ymin = 2^lwr,
                  ymax = 400),
              fill = "green",
              alpha = 0.2,
              inherit.aes = FALSE) + 
  geom_line(data = pi_20_80,
            aes(x = 2^l2pb,
                y = 2^fit),
            linetype = "dashed",
            inherit.aes = FALSE) +
  scale_x_continuous("Bleeding per 10^4 patients",
                     limits = c(50, 400),
                     expand = c(0, 0)) +
  scale_y_continuous("Stroke per 10^4 patients",
                     limits = c(0, 400),
                     expand = c(0, 0)) +
  ggtitle("Fear of Bleeding: 80/100\nGREEN = MEDICATE") +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5))
print(p1)

4.6 Plot Fear Bleeding 80 on 0-100 Original Scale

p2 <- ggplot(tt,
             aes(x = pb,
                 y = ps,
                 fill = lab)) +
  geom_errorbar(aes(ymin = 10^4*(prob.stroke - 1.96*se.stroke),
                    ymax = 10^4*(prob.stroke + 1.96*se.stroke)),
                size = 0.5) +
  geom_errorbar(aes(xmin = 10^4*(prob.bleed - 1.96*se.bleed),
                    xmax = 10^4*(prob.bleed + 1.96*se.bleed)),
                size = 0.5) +
  geom_point(shape = 21,
             aes(size = sqrt(N))) +
  geom_ribbon(data = pi_20_80,
              aes(x = 2^l2pb,
                  ymin = 2^upr,
                  ymax = 400),
              fill = "green",
              alpha = 0.2,
              inherit.aes = FALSE) + 
  geom_line(data = pi_20_80,
            aes(x = 2^l2pb,
                y = 2^fit),
            linetype = "dashed",
            inherit.aes = FALSE) +
  scale_x_continuous("Bleeding per 10^4 patients",
                     limits = c(50, 400),
                     expand = c(0, 0)) +
  scale_y_continuous("Stroke per 10^4 patients",
                     limits = c(0, 400),
                     expand = c(0, 0)) +
  ggtitle("Fear of Bleeding: 80/100\nGREEN = MEDICATE") +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5))
print(p2)

5 Method 2: Deming Regression

NOTE1: Deming regression requires standard errors for each point estimate. I do not have these estimates on the log scale. Hence, for the demonstration purposes, I’m showing regression models on the original scale.

NOTE2: how to do Deming regression prediction interval?

m1 <- lm(prob.stroke ~ prob.bleed,
         data = tt)
summary(m1)

Call:
lm(formula = prob.stroke ~ prob.bleed, data = tt)

Residuals:
       Min         1Q     Median         3Q        Max 
-0.0060634 -0.0024344 -0.0001482  0.0016504  0.0059096 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.002773   0.002342  -1.184 0.263725    
prob.bleed   0.688676   0.124877   5.515 0.000256 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.003809 on 10 degrees of freedom
Multiple R-squared:  0.7526,    Adjusted R-squared:  0.7278 
F-statistic: 30.41 on 1 and 10 DF,  p-value: 0.0002564

m2 <- deming(prob.stroke ~ prob.bleed,
             ystd = se.stroke,
             xstd = se.bleed,
             data = tt)
m2

Call:
deming(formula = prob.stroke ~ prob.bleed, data = tt, xstd = se.bleed,     ystd = se.stroke)

n= 12
                  Coef    se(coef)   lower 0.95   upper 0.95
Intercept -0.003480583 0.001783199 -0.006975588 1.442144e-05
Slope      0.726042116 0.218915292  0.296976028 1.155108e+00

   Scale= 2.014618 

p3 <- ggplot(tt,
             aes(x = prob.bleed,
                 y = prob.stroke,
                 fill = lab)) +
  geom_errorbar(aes(ymin = prob.stroke - 1.96*se.stroke,
                    ymax = prob.stroke + 1.96*se.stroke),
                size = 0.5) +
  geom_errorbar(aes(xmin = prob.bleed - 1.96*se.bleed,
                    xmax = prob.bleed + 1.96*se.bleed),
                size = 0.5) +
  geom_point(shape = 21,
             aes(size = sqrt(N))) +
  geom_abline(intercept = m1$coefficients[1],
              slope = m1$coefficients[2]) +
  geom_abline(intercept = m2$coefficients[1],
              slope = m2$coefficients[2],
              linetype = "dashed") +
  ggtitle("Least Squares (solid line) vs.\n Deming (dashed line) Regression") +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5))
print(p3)

6 Session Information

sessionInfo()

R version 4.0.3 (2020-10-10)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 17763)

Matrix products: default

locale:
[1] LC_COLLATE=English_United States.1252  LC_CTYPE=English_United States.1252   
[3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C                          
[5] LC_TIME=English_United States.1252    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] MethComp_1.30.0   deming_1.4        plotly_4.9.3      ggplot2_3.3.3    
[5] data.table_1.13.6

loaded via a namespace (and not attached):
 [1] pillar_1.4.7      compiler_4.0.3    tools_4.0.3       boot_1.3-25      
 [5] digest_0.6.27     evaluate_0.14     jsonlite_1.7.2    lifecycle_0.2.0  
 [9] tibble_3.0.5      gtable_0.3.0      nlme_3.1-151      viridisLite_0.3.0
[13] lattice_0.20-41   pkgconfig_2.0.3   rlang_0.4.10      cli_2.2.0        
[17] rstudioapi_0.13   crosstalk_1.1.1   yaml_2.2.1        xfun_0.20        
[21] coda_0.19-4       withr_2.4.0       dplyr_1.0.3       httr_1.4.2       
[25] knitr_1.30        generics_0.1.0    vctrs_0.3.6       htmlwidgets_1.5.3
[29] grid_4.0.3        tidyselect_1.1.0  glue_1.4.2        R6_2.5.0         
[33] fansi_0.4.2       rmarkdown_2.6     farver_2.0.3      purrr_0.3.4      
[37] tidyr_1.1.2       magrittr_2.0.1    scales_1.1.1      ellipsis_0.3.1   
[41] htmltools_0.5.0   assertthat_0.2.1  colorspace_2.0-0  labeling_0.4.2   
[45] lazyeval_0.2.2    munsell_0.5.0     crayon_1.3.4